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ABSTRACT 



We present a near-infrared extinction study of the dark globule B335, a pro- 
tostellar collapse candidate, using data from HST/NICMOS and the W.M. Keck 
Observatory. These data allow a new quantitative test of the "inside-out" collapse 
model previously proposed to explain molecular line profiles observed toward this 
region. We find that the shape of the density profile is well matched by the col- 
lapse model, but that the amount of extinction corresponds to larger column 
densities than predicted. An unstable Bonnor-Ebert sphere with dimensionless 
outer radius £ max = 12.5 ±2.6 provides an equally good description of the density 
profile, and is indistinguishable from the collapse model over the range in radius 
sampled by the extinction data. The bipolar outflow driven by the embedded 
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young stellar object has an important effect on the extinction through the core, 
and modeling the outflow as a hollowed-out bipolar cone of constant opening 
angle provides a good match to the observations. The complete exinction map 
is well reproduced by a model that includes both infall and outflow, and an ad- 
ditional 20% dispersion that likely results from residual turbulent motions. This 
fitted model has an infall radius of = 26 ± 3" (0.031 pc for 250 pc distance), 
and an outflow cone semi-opening angle of a = 41 ±2°. The fitted infall radius is 
consistent with those derived from molecular line observations and supports the 
inside-out collapse interpretation of the density structure. The fitted opening an- 
gle for the outflow is slightly larger than observed in high velocity CO emission, 
perhaps because the full extent of the outflow cone in CO becomes confused with 
ambient core emission at low velocities. 

Subject headings: ISM: globules — ISM: individual(B335) — dust, extinction - 
ISM: jets and outflows — stars: formation 



1. Introduction 

The near-infrared color excess of extincted background stars can reveal the dust dis- 
tribution in a dense core and provide information on physical structure. Historically, this 
technique has been hampered in small, dense regions by limited instrumental sensitivity 
(Jones et al. 1980, 1984), but the development of large format array cameras has sparked 
renewed interest (Lada et al. 1994, Alves et al. 1998, Lada, Alves & Lada 1999, Alves, Lada 
& Lada 1999, 2001). The unprecedented sensitivity available with the Near Infrared Cam- 
era and Multi Object Spectrometer (NICMOS) instrument on the Hubble Space Telescope 
(HST) opens up a new regime. In particular, the high extinctions of dense cores that form 
the immediate environs of candidate low mass protostars are now accessible. 

Molecular line surveys of nearby dark clouds have identified a large number of low 
mass dense cores (Myers, Linke & Benson 1983) of which roughly half are associated with 
young stellar objects detected by IRAS (Beichman et al. 1986). These observations have 
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contributed to the general paradigm that at least some low-mass stars like the Sun form by 
the gravitational collapse of isolated, dense cores. In the standard theory, originating with 
Shu (1977), the dynamical collapse to a star/disk system occurs from the inside-out. First, 
an n oc r~ 2 distribution is established prior to collapse as the cloud core loses magnetic 
and turbulent support through ambipolar diffusion and relaxes to a balance between gravity 
and thermal pressure. Collapse is initiated at the center where the density is highest, and 
a wave of infall propagates outward at the sound speed. Conditions inside the infall radius 
asymptotically approach free-fall, with n oc r~ L5 and v oc r -0 - 5 . Modifications to the simple 
theory have been made to account for more realistic initial conditions, including slow rotation 
(Tereby, Shu & Cassen 1984) and residual magnetic fields (Galli & Shu 1993, Li & Shu 1997). 
A similar radial density profile is also produced by the competing self-similar collapse theory 
of Larson (1969) and Penston (1969). However, the Larson-Penston flow predicts far larger 
infall velocities than the Shu theory, and has proved unsuccessful at explaining molecular 
line observations (see e.g. Zhou 1992). 

Unambigous identification of inside-out collapse motions has proved difficult for several 
reasons. First, the fact that collapse initiates in the center of a core means that it can be 
detected only with observations of sufficiently high angular resolution. Further difficulty 
comes from the fact that the velocities associated with infall motions are small, comparable 
to intrinsic line widths and much smaller than the velocities of the bipolar outflows powered 
by young stellar objects. The presence of localized, redshifted self- absorption in spectral 
lines of moderate optical depth provides one signature of infall motions, but this signature is 
not unique. Scenarios other than collapse can produce the same spectral line profiles. The 
most promising claims for collapse have been inferred from redshifted self-absorption, and 
consequently these claims are often controversial (see the review of Myers, Evans & Ohashi 
2000). 

The dense core in the B335 region is generally recognized as the best protostellar collapse 
candidate. The B335 region contains an isolated, roughly spherical globule at a distance of 
roughly 250 pc (Tomita, Saito & Ohtani 1979). Deeply embedded in the globule is a young 
stellar object discovered at far-infrared wavelengths by Keene et al. (1983) and detected only 
at A > 60 /iin by IRAS. Observations of high velocity molecular line emission have shown 
that the young stellar object drives a bipolar outflow with semi-opening angle of ~ 25 ± 5° 
(Cabrit, Goldsmith & Snell 1988, Hirano et al. 1988, 1992). Detailed radiative transfer 
models based on the theory of inside-out collapse provide good fits to many molecular line 
profiles of dense gas tracers observed at 10"-30" resolution (Zhou et al. 1993, Choi et al. 
1995). These models imply an infall radius of 0.03 pc (25"), an age of 1.3 x 10 5 yrs, and 
a central mass of 0.37 M Q . While the agreement of the model with few parameters and 
the observed molecular line profiles strongly support inside-out collapse, the model is not 
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unique. For example, in one extreme scenario discussed by Zhou et al. (1993), an unrelated 
foreground cloud of just the right excitation properties absorbs emission from a static core 
and exactly reproduces the infall signatures. Recent molecular line observations of B335 
made with the higher angular resolution of interferometry are also starting to show perhaps 
important discrepancies with the predictions of inside-out collapse (Wilner et al. 2000). 

An extinction map of sufficient depth can probe the density structure of B335 and test 
the claim for inside-out collapse in a way that is free from many of the problems that plague 
molecular line studies, including excitation effects and spatially varying abundances of the 
molecular tracers. B335 is fortuitously placed for an extinction study, located at Galactic 
coordinates I = 44.9, b = —6.6, which is a direction with a very high density of background 
stars with uniform properties. Previous optical extinction studies of B335 were not sensitive 
to the structure of the dense core anywhere close to where the molecular line data suggest 
collapse (see Bok & McCarthy 1974, Dickman 1978, Tomita et al. 1979 and Frerking et al. 
1987). In this paper, we present a near-infrared extinction study of B335 using data from 
HST/NICMOS, supplemented by data from the W.M. Keck Observatory. Section 2 briefly 
reviews the technique and describes the observations and data reduction. Section 3 presents 
an analysis of the extinction data in the context of several physical models for the B335 
density structure, and comparisons with the results of previous studies. Simulations of the 
observations are also presented, using a Monte Carlo approach detailed in the appendix. 
Section 4 summarizes the main conclusions. 



2. Observations and Data Reduction 

2.1. Review of the Near-Infrared Color Excess Technique 

Measurements of the extinction of stars background to a dense core provides estimates 
of column density along many pencil beams. These estimates can be used to determine 
the overall structure of the obscuring core. The extinction measurements are best done 
at infrared wavelengths, where absorption due to dust is much less than in the optical, 
thereby probing regions of higher (visual) extinction. The basic method is to measure the 
near-infrared color excess for each background star: 

E(H - K) = (H - K) obscrvcd — (H — K), , (1) 

where (H — K)* is the intrinsic color of the star. The color excess represents the differential 
optical depth between the two near-infrared filters, which is directly proportional to dust 
column density. The observed color excess may be converted to an equivalent visual extinc- 
tion using the standard reddening law (e.g. Rieke & Lebosfsky 1985), and then converted to 
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gas column density by adopting a gas-to-dust ratio (e.g. Bohlin et al. 1978). There are some 
uncertainties in both of these conversions. The gas-to-dust ratio in protostellar cores is the 
subject of some debate; at low extinction there is some evidence to suggest that it is constant 
along different lines of sight in the galaxy (Bohlin et al. 1978), but at high extinction it is 
unclear whether the gas-to-dust ratio remains constant. For example, if grain growth occurs 
in the high extinction regions as a result of molecular depletion and agglomeration, then one 
might expect lower extinction for a given column density. The reddening law is probably less 
problematic and nearly universal at near-infrared wavelengths, given the small grain sizes 
present in protostellar envelopes (Mathis 1990). 

In practice, the color excess method does not require a knowledge of the intrinsic colors 
of the individual extincted stars. A statistical correction to their observed colors may be ob- 
tained empirically from the background stellar population, provided it is sufficiently spatially 
uniform. For stars of spectral types AO to late M, in the main sequence or the giant loci, the 
intrinsic colors lie in the narrow range < (H — K)* < 0.3 (Koornneef 1983). Therefore, at 
high extinctions, little error is introduced by simply assuming a mean background star color 
derived from nearby control fields. 

2.2. NICMOS Observations 

The B335 globule was imaged with the NIC3 camera on HST (Cycle 7— NICMOS, 
GO-7843, 1998 June 23) using two broadband filters, F160W and F222M, similar to the 
usual H band and K band filters that are matched to the 1.6 fim and 2.2 fim atmospheric 
windows. The NIC3 camera has 256 x 256 pixels with size C/20 giving a 51'/ 2 field of view. 
The observations were obtained during the second of two periods of a few weeks when the 
HST secondary mirror was adjusted to allow the NIC3 camera to be in focus (the "NIC3 
Campaigns" ) . 

The observing program had two parts: (1) a deep 3x3 mosaic of the central region of 
high extinction, taking 6 orbits, and (2) a radial strip of images extending about 4' from the 
center, taking 1 orbit, to reach the background stellar population and provide a sample of 
stars to transform NICMOS magnitudes to a standard photometric system. The length of the 
strip was the maximum allowed before suffering the heavy overhead of guide star reaquisition. 
(The strip included the center of B335 as insurance in the event that the deeper observations 
were unsuccessful.) Figure 1 shows the observed fields overlayed on a wide field optical view 
of the region from the Digital Sky Survey. We note that the NIC3 orientation angle on the 
sky and stepping direction of the radial strip were given by the spacecraft orientation on 
the day of observation. Advance specification of a particular orientation would have created 
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an untenable scheduling constraint; the orientation obtained was based solely on where the 
spacecraft happened to be during the period when the NIC3 camera was in focus and B335 
was accessible. 

The central region was observed with a standard nine-point spiral dither pattern with 
10" offsets, and a total exposure time of 896 seconds per dither. The resulting square mosaic 
covers 72". The strip images were taken in one-chop mode, an exposure time of 56 seconds 
per chop, with an offset of 3" between chops, and roughly 3"7 overlap between successive 
images. This resulted in a strip of coverage 55" wide that extends 247" from the mosaic 
center. The data were taken in the NICMOS multiaccum mode, for which an initial readout 
of the detector array is followed by non-destructive readouts during the course of a single 
integration. The intermediate readouts can be used to identify and to remove the effects of 
cosmic ray hits and saturated pixels. 

The NICMOS data were obtained from the HST archive in FITS format both unpro- 
cessed and already processed through the STScI pipeline (CALNIC A & B). Initial inspection 
showed the STScI pipeline did not adequately deal with a variety of instrumental effects, and 
so we reprocessed the data in IRAF using the NICREDUCE software developed by Brian 
McLeod (McLeod 1997, Lehar et al. 1999). The NICREDUCE package provides superior 
and more versatile handling of sky and dark current subtraction, as well as cosmic ray hits 
and other defects. 

The bottom rows of the NICMOS images suffered from a small amount of vignetting. 
The effect was most pronounced in the F222M images, where about 2 / .'4 were unusable. The 
overlap of images in the dither pattern and in successive images of the strip was sufficient to 
ensure that the vignetting resulted in only a small loss in sensitivity in the affected regions, 
and no gaps in the spatial coverage. 

The dithered images of the central region of B335 suffered an additional complication 
due to the presence of a very bright star (H magnitude of 9.0) on the west side of the 
mosaic center. The enormous count rate from the star caused the image reduction software 
to produce a small depression in the continuum level of the quadrant of each image that 
contained the star. This deviation was fixed by adding a constant to the relevant quadrant 
so that for each image the median pixel value was uniform in all four. Ghost impressions of 
the star in successive exposures were interactively masked in the affected images and were 
not included in the final mosaic. 
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2.2.1. Photometry 

Stars were identified by visual inspection of the images, as several attempts at auto- 
matic star finding were compromised by the effects of diffraction spikes and rings around 
the brighter stars. Photometry on each star was performed with the suite of routines in 
the apphot package. A series of apertures from Of! 2 to 0'/5 in radius was used for each star. 
Aperture corrections were then performed using a minimum of four isolated and bright stars 
from the relevant image. Each aperture was used to calculate a magnitude and associated 
error as a check to identify any unusual behavior. The measured count rates were corrected 
to a 1" radius aperture, and then multiplied by 1.075 as detailed in the NICMOS Photo- 
metric Calibration CookBook (available from STScl). The instrumental magnitudes were 
then transformed to the Vega scale using constants for each filter provided in the HST Data 
Handbook (Version 4, Dec 1999). For stars that were observed in a region of overlap between 
strip images, weighted mean magnitudes were calculated. For stars observed in the overlap 
region of the central mosaic and the first strip image, the value from the more sensitive 
central mosaic was taken. The complete NICMOS photometry list contains 1026 stars de- 
tected with the F160W filter and 471 stars detected with the F222M filter with magnitude 
error < 0.1. The NICMOS limiting magnitudes are spatially dependent and are depicted in 
Figure 2. In the central 32", with the most integration time, the limiting magnitudes are 
H=24.8 and K=21.5 (in the CIT system). 



2.3. Keck Observations 

Observations of B335 were made on 1998 July 10 using NIRC (Matthews & Soifer 1994) 
on the Keck I telescope. These observations were made for two reasons: (1) to obtain a 
better study of the stellar background population around B335, and (2) to transform the 
NICMOS F160W and F222M magnitudes to a standard system. The NIRC has 256 x 256 
pixels with size Of! 15, giving a 38'.'4 field of view. Each basic observation consisted of nine 
dithered 5 second integrations, in a 3 x 3 pattern with 3" offsets, which were mosaiced 
together. Observations were made with J, H, and K filters. The program had three parts: 
(1) nine images of the central region of B335, approximately mimicking the sub-images of 
the NICMOS central mosaic, (2) five images leading away from the center of B335, at each 
position of a NICMOS strip image, and (3) four additional fields 10' away from B335, to 
study the characteristics of the stellar background population. These regions are marked on 
Figure 1. 

The data reduction and calibration was done in IRAF following the standard dark/flat /sky 
subtraction procedures, using a running flat. This procedure works well when the background 
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emission does not change rapidly. For some periods during the observations, the sky emission 
in the H band did change relatively quickly, and the background subtraction for observations 
during these periods is not perfect. The typical seeing during the night was Cf/5, though there 
were periods when it was better and also much worse. 

2.3.1. Photometry 

The majority of stars in the Keck images of B335 could be identified with reference to 
the NICMOS data, since the Keck images are all less deep than the NICMOS equivalents 
and the field of view slightly smaller. For the Keck fields not observed with NICMOS, all 
stars were again identified by visual inspection. 

Photometry was performed for all non-saturated stars using the apphot package. Again, 
a series of increasing radial apertures was used for each star, and each aperture was used 
to calculate a magnitude. This approach allowed a confidence level in the measurement to 
be assigned to each star, based on whether the distribution of measurements from different 
apertures was consistent with the predicted uncertainty. If the spread of the various magni- 
tudes was greater than 0.2 magnitudes, then the measurement was identified as low quality 
(18% of the stars detected at H and 17% of the stars detected at K). These low quality stars 
have point-spread-functions that differ from the norm, mainly from those periods of poor 
and varying seeing conditions. The instrumental magnitudes were transformed to the CIT 
system through observations of the HST infrared standards SJ9164 and SJ9182 (Persson et 
al. 1998) and the UKIRT standard FS35 (Casali & Hawarden 1992). The internal disper- 
sion in the zero points measured for these stars was 0.02 magnitudes, and no atmospheric 
extinction correction could be discerned. The Keck photometry list contains a total of 297 
stars detected at H, and 340 stars detected at K, in the overlap region between the Keck 
and NICMOS datasets. 

The four Keck fields 10' away from B335 contain typically 100 detected stars per square 
arcminute with H< 20. The four fields display uniform properties and no discernable red- 
dening; the star counts in each image are identical to within Poisson errors, and the mean 
(H — K) color of the fields are all in the range 0.07 to 0.11, consistent with the intrinsic 
colors of giant and main sequence stars. A detailed study of the background stellar pop- 
ulation is presented in Section 3.1. The limiting magnitudes of the Keck observations are 
approximately H=19.7 and K=19.6. 
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2.4. Transformation of NICMOS Magnitudes 

For extinction calculations, the NICMOS F160W and F222M magnitudes are trans- 
formed to standard H and K magnitudes. Linear least squares fits (not including the lower 
quality Keck data) give slopes of -0.05 for (H-F160W) vs. F160W, -0.03 for (K-F222M) 
vs. F222M, and 0.99 for (H-K) vs. (F160W-F222M), where the formal uncertainties in the 
fitted slopes are ~ 0.01. The formal uncertainties are underestimates of the true uncertain- 
ties since the magnitude comparison slopes are strongly biased by the few brightest stars, 
and the color comparison slope is biased by the few reddest stars. Experimentation with 
subsets of the stars resulted in variations in the fitted slopes at the 0.05 level. We conclude 
that these data provide no evidence for non-zero slopes or color terms in the transformation. 
We calculated constant offsets for the magnitude conversions for each filter, using a weighted 
mean of (H-F160W) and (K-F222M) values. The result is H = F160W+ (0.035 ±0.003), 
and K = F222M + (0.052 ± 0.003), where the uncertainties represent the standard error in 
the mean. Figure 3 shows the scatter plot of (H-K) vs. (F160W- F222M) with the derived 
color transformation superimposed. For the (H — K) colors of interest, in the range to 3, 
the differences between the two photometric systems are remarkably small. 

Support for the accuracy of these transformations comes from published ground-based 
photometry of the one very bright star observed by Krugel et al. (1983), who measured 
Johnson H and K fluxes of 310 ±25 mJy and 605 ±45 mJy, which correspond to magnitudes 
of H = 8.89 ± 0.07 and K = 7.58 ± 0.07, and color (H - K) = 1.31 ±0.11. These values 
compare very favorably with the transformed NICMOS magnitudes of H = 9.039 ± 0.001 
and K = 7.675 ±0.001, and color (H — K) = 1.364 ±0.001. Note that systematic errors limit 
the accuracy of the NIC3 photometry to a few percent, and are not reflected in the random 
uncertainties quoted here. 

3. Results and Analysis 

3.1. The Background Population Mean Color and Dispersion 

To check the uniformity of the background stellar population, the starcounts and mean 
colors of the four Keck off fields were compared with those for the two NICMOS strip images 
furthest away from the center of B335. The colors in the different images are consistent 
within the observed dispersions. The fifth NICMOS strip field has a mean color and standard 
deviation of H — K = 0.17, <x(H — K) = 0.18, very similar to that found for the combined 
Keck off-fields: H-K = 0.10, a(R - K) = 0.13. The fourth NICMOS strip field shows a 
slightly redder mean color (H — K = 0.26) and higher dispersion (cr(H — K) = 0.23) than 
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the other fields; most likely, this field is not far enough away from the center of B335 to 
be completely free of its extinction, and this field was excluded from our characterization of 
the background stellar population. Figure 4 presents a histogram of the (H — K) colors of 
the "background" stars in 0.1 magnitude bins (excluding the Keck stars with lower quality 
photometry). The (H — K) mean color is 0.13 magnitudes, with standard deviation 0.16 
magnitudes. Also shown for comparison is a Gaussian distribution with the same mean 
and standard deviation (normalized so that the area under each curve is equal). A color- 
magnitude diagram constructed for this sample shows no sign of any important systematic 
trends, e.g. there is no tendency for fainter stars to be redder. The starcounts in the 
fifth NICMOS strip field are indistinguishable from those in the Keck off-fields, taking into 
account the different different size fields of view. 



3.2. NICMOS Central Mosaic Image 

Figures 5 & 6 show the NICMOS F160W and F222M images of the B335 region. The 
central mosaic images in Figure 5 may be compared to H band and K band images of a 
smaller part of the central region of B335 obtained over three nights at Keck by Hodapp 
(1998) using NIRC. In the region of overlap, the F222M image shows the same stars as the 
K band image, but the overall sensitivity of the NICMOS image is better and the point- 
spread-function much sharper, leading to better separation of nearby stars. The F160W 
image shows many more stars than the corresponding H band image. In contrast with the 
ground-based data, essentially all stars detected with the F222M filter are also detected with 
the F160W filter because of the low background and greater sensitivity of the space-borne 
NICMOS at the shorter wavelength. 

Figure 7 conveys much of the information from the central mosaic of NICMOS obser- 
vations. The left plot is a pseudo-image in that the two axes are spatial, the apparent 
brightness of a star determines its "size" , and the value of its (H — K) color determines its 
"color". The right-hand plot displays the radial dependence of the (H — K) colors out to 180" 
from the center. The origin in both plots is given by the position of the embedded young 
stellar object, as determined by interferometric observations of millimeter continuum emis- 
sion (Wilner et al. 2000); the uncertainty in the registration of the NICMOS images and the 
millimeter peak is estimated to be about 1". Based on an inspection of the extinctions, the 
stars in the right hand plot have been marked according to their spatial location: those that 
lie in the North or South quadrants are marked in red, and those that lie in the East or West 
quadrants are marked in blue (the quadrant boundaries are marked in the pseudo-image as 
dashed lines). Defining 9 to be the angle subtended between the vector to a star and the axis 
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of the outflow (which runs very close to east-west), this division may be expressed simply 
as 9 > 45° or 9 < 45°. The NICMOS central mosaic image contains 239 stars detected at 
H, 121 stars detected at K, and 119 stars that are detected in both filters. The detected 
star closest to the embedded young stellar object is at 14" radius. The reddest detected star 
is at 17" radius and has an (H — K) color of 4.0, corresponding to over 60 magnitudes of 
equivalent visual extinction. 

The major features of Figure 7 are: (1) a strong gradient in the (H — K) colors as 
one approaches the protostar; (2) fewer stars are detected close to the center, despite bet- 
ter sensitivity; (3) reddening of the stars turns off into background noise beyond the 125" 
(0.15 pc) outer radius of dense gas determined from molecular line data (Zhou et al. 1990); 
(4) there is a marked asymmetry in the stellar colors associated with the bipolar outflow- 
stars viewed within ~ 45° of the outflow axis have systematically lower colors; (5) the dis- 
persion of the (H — K) colors at a given radius and in the quadrants away from the outflow 
is ~ 20% x E(H — K) greater than that of the background population. 

3.3. Theoretical Models of the B335 Density Distribution 

We describe several theoretical models for the density distribution of B335 and evaluate 
the success of these models in light of the extinction data. The models considered are not 
meant to comprise an exhaustive list. They include the inside-out collapse models previously 
proposed for B335, Bonnor-Ebert spheres, and simple power law descriptions. We also discuss 
briefly the effects of the bipolar outflow and turbulent motions in B335. 

3.3.1. Inside- Out Collapse Models 

There is a substantial literature of theoretical work by Frank Shu and colleagues describ- 
ing the density structure of collapsing dense cloud cores, starting from the initial condition of 
a singular isothermal sphere. Taking parameters appropriate for the B335 region, we briefly 
describe the pre-collapse state, the spherical inside-out collapse, and the modifications to 
spherical collapse expected for slow initial rotation or weak magnetic fields. 

Pre-collapse Configuration: The initial state is that of a singular isothermal sphere; the 
density falls off as p oc r~ 2 , where r is the radius, with the normalization of the density 
determined by the effective sound speed. For B335, Zhou et al. (1990) derive a value for the 
effective sound speed of a = 0.23 km s -1 from molecular line data, which corresponds to a 
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kinetic temperature of 13 K. The static initial condition is therefore given by: 

a 2 

Pstatic(r) = (2) 

= 1.33 x 10- 20 gem" 3 ( £ Y(^—) 2 ( 3 ) 

n^ 2 (r) = 2.80 x 10 3 cm" 3 ( ^ V ( —!—] (4) 

where the conversion to a molecular hydrogen number density assumes 0.7 for the the mass 
fraction in hydrogen. 

Spherically Symmetric Collapse: The simplest scenario is the spherically symmetric 
collapse (Shu 1977). In this model, a spherical wave of collapse propagates outwards from 
the center at the effective sound speed. The radial distance that the wave has traveled, 
sometimes called the infall radius, is the only additional parameter that defines the density 
distribution. Inside the infall radius, conditions approach free fall, with the density taking 
the asymptotic form p oc r -3 / 2 . By reproducing essential features of several molecular line 
profiles with this model, Zhou et al. (1993) claimed evidence for inside-out collapse in B335. 
Choi et al. (1995) made more detailed radiative transfer models to match the Zhou et al. 
(1993) observations of H 2 CO and CS molecules, determining the infall radius to be 0.030 pc 
(25"). Note that this derived infall radius falls entirely within the deep central mosaic of 
NICMOS coverage. 

Rotation: The effect of slow rotation on the collapse has been solved using a perturbative 
approach from the spherically symmetric case (Terebey, Shu & Cassen 1984). The result is 
a distortion to the shape of the collapse wave of order 4 x 10~ 4 r 2 , where r is the number 
of times that the cloud has rotated since collapse initiated (this is also equivalent to the 
ratio of the infall radius to the centrifugal radius). For B335, Frerking et al. (1987) found 
an angular velocity of 1.4 x 10 -14 s _1 through C 18 observations. This indicates a value of 
r 2 ~ 3 x 10~ 3 , and a distortion in the wave by only ~ 10~ 6 — an entirely negligible effect. 
The density distribution is more strongly distorted, with the effect largest near the center 
of collapse. For radii greater than a tenth of the infall radius (i.e. r > 2'/5) — already well 
within the innermost radius where stars are detected — the correction is always less than 
of order r 2 , which is also completely negligible. Such slow rotation also has an unimportant 
effect on the model line profiles towards the center of B335 (Zhou 1995). We therefore do 
not include rotation in our models. 

Magnetic Fields: The effect of a weak magnetic field on the inside-out collapse of a 
cloud has been studied by Galli & Shu (1993) and Li & Shu (1997). The collapse wavefront 
and the density distribution are again distorted, in a similar way to the rotational model. In 
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the magnetic case, the parameter r is given by the ratio of the infall radius to the magnetic 
radius: 

Galli & Shu (1993) showed that the field effects only a small perturbation to the spherically 
symmetric solution outside a radius of r > 0.2 r 4 / 3 R mi . Figure 7 shows that variation in 
color between the two spatial populations of stars can be as large as a factor of two, well 
beyond the infall radius of 25" derived by Choi et al. (1995). Since this certainly cannot 
be construed as a "small perturbation" at these radii, we can place a lower limit on the 
magnetic field that would be needed to explain the asymmetry in the context of this model. 
Taking r/i? in f ~ 1, we find that: 

B «^ 50 (a») 2 (^)"" G < 7 > 

Such a minimum magnetic field strength is implausibly large. For example, Crutcher & 
Troland (2000) measured a line-of-sight magnetic field of B ~ 11 p,G in the starless pro- 
tostellar core L1544 using the Zeeman effect in lines of OH, and this is much higher than 
the typical upper limits found for other dark clouds (Crutcher et al. 1993). In this context, 
it therefore seems very unlikely that a magnetically modified inside-out collapse could be 
responsible for the large asymmetry observed in the B335 extinction data. 



3.3.2. Bonnor-Ebert Models 

Bonnor-Ebert models are pressure confined isothermal spheres, for which the solution 
remains physical at the origin (Ebert 1955, Bonnor 1956). In common with the singular 
isothermal sphere, the initial condition for inside-out collapse, they are solutions of a modified 
Lane-Emden equation (Chandrasekhar 1967): 

^(«*)-«pH». (s, 

where £ = (r/R ) = (r/a)y/ATvGp c is the dimensionless radius, and = — \n(p/p c ) is a 
logarithmic density contrast, with p c the (finite) central density. Unlike the singular solution, 
the Bonnor-Ebert solutions do not diverge at the origin. Instead, the boundary conditions 
are that the function ip and its first derivative are zero at the origin. 
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The above equation can be solved by division into two first order equations which can 
then be tackled simultaneously using numerical techniques (in our case a 4th order Runge- 
Cutta method). There is a family of solutions characterized by a single parameter — the 
dimensionless outer radius of the sphere, £ max - For a given sound speed and a particular 
choice of the shape of the density curve (i.e. £ m ax), there is one additional degree of freedom, 
the physical scale of the model. This additional degree of freedom may be removed by 
implementing an additional constraint, for example by fixing the outer radius, the total 
mass, or the external pressure. 

Configurations with dimensionless outer radius £ max > 6.5 are unstable to gravitational 
collapse (Bonnor 1956). The gravitational collapse of Bonnor-Ebert spheres has been studied 
numerically by Foster & Chevalier (1993). They find that the flow asymptotically approaches 
the Larson-Penston solution at the origin at the time of and prior to the formation of a central 
core, but emphasize that these large early accretion rates last for a very short time. If the 
cloud is initially very centrally condensed (i.e. £ max 3> 6.5) the later stages of infall closely 
resemble Shu's inside-out collapse of a singular isothermal sphere. 

Recently, Alves, Lada & Lada (2001) used the color excess method to study B68, a 
starless core, with data from the ESO VLT. The B68 system is much less centrally condensed 
than B335, and has much lower column density at the center (only about 30 magnitudes of 
equivalent visual extinction). The Alves et al. (2001) extinction map of B68 suggests that the 
density structure is well described by a Bonnor-Ebert sphere with dimensionless outer radius 
slightly in excess of critical: £, m ax — 6.9 ± 0.2. Although the B335 core is at a much later 
evolutionary stage than B68, with an embedded protostar, we consider the Bonnor-Ebert 
sphere since it describes B68 so successfully. 

3.3.3. Bipolar Outflow and Turbulence 

The outflows from young stellar objects can have a substantial impact on their surround- 
ings. A bipolar outflow in B335 was discovered in CO by Frerking & Langer (1982). Later 
studies in various tracers by Cabrit, Goldsmith & Snell (1988), Hirano et al. (1988, 1992), 
Chandler & Sargent (1993), and Wilner et al. (2000) show that on arcsecond to arcminute 
scales, the outflow is well collimated with a semi-opening angle of about 25 ±5°, and with the 
axis lying very nearly in the plane of the sky. The outflow orientation is close to east-west, 
with position angle roughly 100° in the coordinate system of the NICMOS images. Inside 
the region affected by the outflow, the density structure is modified as the core material is 
swept out by the flow, but the exact form of the density field associated with the outflow is 
not known. 
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The bipolar outflow may be an important driver of turbulent motions. Zhou et al. (1990) 
determined a turbulent component to the effective sound speed in B335 of 0.08 km s -1 . Cou- 
pling of the velocity and density fields via the continuity equation suggests that the pertur- 
bations to the velocity field should induce density contrast perturbations, i.e. 5 ~ [(1/3) < 
f/a> 2 ] 1 ^ 2 ~ 0.2. The resulting perturbations to the column densities (and extinctions) will 
depend on the size scale of the velocity and density field perturbations. The perturbations 
to the column density will be largest when the size scale is a factor of a few smaller than the 
dense core. The extinction map is consistent with this possibility since the observed colors 
of neighboring stars are generally very similar, yet there exists a ~ 20% dispersion for stars 
that are not neighbors but are located at similar radial distances from the protostar (see 
Figure 7). The effect of the column density perturbations are inevitably reduced by inte- 
gration along the line-of-sight, but the observed residual turbulence most likely contributes 
significantly to the dispersion in the observed colors. 



3.4. Fitting Model Parameters and Evaluating Fit Quality 

We evaluate the goodness of fit for these various model density distributions by calcu- 
lating a reduced x 2 , defined as: 

2 J_ V (^ E(H-K)f IC -E(H-K)r dcl \ 2 , Q . 

" N-m -NIC ' W 



N — m ^-^ \ a % 

where m is the number of free parameters in the model being fitted, and the sum extends over 
a particular subset of N stars taken from the total number detected in both NICMOS filters. 
The values of the observed color excess E(H — K) 4 NIC and the uncertainty af 10 are calculated 
using the properties of the unreddened stellar background population (see Section 3.1, and 
Figure 4), assuming each star to have an intrinsic color of (H — K) B q = 0.13, with an 
uncertainty of <tbg = 0.16: 

E(H-K)f IC = (H-K)f IC -(H^K) BG 

= (H — K)j NIC — 0.13 , (10) 



^NIC I 2 i _2 

i = \ a i + CT BG 



= ^a 2 + 0.16 2 , (11) 

where <7j is the uncertainty in the observed (H — K) color of a given star. This procedure 
implicitly assumes that the photometric errors (<7j) and the intrinsic (H — K) colors both 
have Gaussian distributions. Though the latter condition is not strictly true, Monte Carlo 
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trials show that the xl values calculated with the above recipe agree very closely with those 
produced when the intrinsic color of each star is instead chosen randomly from the Gaussian 
distribution with the same mean and standard deviation as the sample of background colors. 
We conclude that the results of the fitting are not sensitive to this approximation. 

Since the models are in general non-linear in the fitting parameters, we analyze the 
uncertainty in the best-fit model parameters using a Monte Carlo technique known as the 
Bootstrap method (Press et al. 1992). To describe the method, let us represent the dataset 
that is being used in the fitting by S. The first step is to construct a new dataset of 
equal size S', by randomly selecting N times from the original dataset. This new sample 
is then analyzed in the same way as the original dataset, and the fitted model parameters 
recorded. This process is then repeated n times, where n is sufficiently large that the 
resulting distribution of best-fit model parameters is insensitive to its exact value (for our 
models, typically n = 200). The standard deviation of this distribution provides an estimate 
for the uncertainty of the parameters that best fit the original dataset. Table 1 summarizes 
the results from the various \ 2 analyses and these results are discussed below. Figure 8 
shows a schematic diagram of the spatial coverage of the NICMOS observations relative to 
the important size scales in B335. 



3.5. Excluding the Outflow (Fits A, B & C) 

The asymmetry in the observed colors appears to be associated with the bipolar outflow. 
Since we do not know from theory the effect of the outflow on the material distribution, we 
initially restrict attention to the quadrants that are located away from the flow, 6 > 45°, 
where is the angle subtended between the vector to a star and the axis of the outflow. To 
avoid noise near the boundary of the globule, we also restrict the fitting region to r < 100". 
In Fit A of Table 1, we use a power law model for the density distribution. In Fit B we use 
the Shu inside-out collapse model, and in Fit C we use a Bonnor-Ebert sphere. 

The best fitting single power law model (Fit A) to the dense core region unaffected by 
the outflow has an index of p — 1.91 ± 0.07. This value is well matched to the r -2 envelope 
of the isothermal sphere that represents the initial condition of the inside-out theory. The 
inside-out collapse solution (Fit B), which is similar to a broken power law with fixed indices, 
provides a somewhat better description than can a single power law model with arbitrary 
index. Remarkably, the fitted infall radius of R in f = 28 ±3" (0.034 ±0.004 pc) agrees within 
uncertainties with those calculated by Zhou et al. (1993) and Choi et al. (1995) by modeling 
molecular line profiles. That the minimum xl value for these fits are near 3 instead of unity 
reflects the ~ 20% x E(H — K) greater dispersion in the observed colors than in the intrinsic 
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colors of the background stars. 

A Bonnor-Ebert sphere provides an equally good description of the density structure in 
B335 (Fit C). Figure 9 shows a comparison of the (H — K) radial profiles predicted by the 
Bonnor-Ebert sphere in Fit C, and the inside-out collapse model in Fit B. The plot shows 
that the two models are practically indistinguishable over the range in radius where there are 
stars to fit (14"-100"). At large radii, the Bonnor-Ebert solution has slope close to a power 
law index of 2.0, and at smaller radii, the slope flattens. In effect, the best fit Bonnor-Ebert 
sphere is one that mimics the inside-out collapse solution over the range in radius for which 
there is data. The models diverge only at such small radii that extinction data of sufficient 
depth to discriminate between them is currently impossible to obtain. One discriminant may 
be sensitive, high resolution observations of dust emission, which is strongest at these small 
radii where extinction data is lacking. Regardless, the value of £ max = 12.5 ±2.6 for the best 
fit Bonnor-Ebert sphere is well in excess of the critical value of 6.5, and it is therefore an 
unstable configuration that should undergo gravitational collapse similar to Shu's inside-out 
solution (Foster & Chevalier 1993). 



3.6. Modeling the Outflow 

The simplest model we imagine to simulate the effects of the bipolar outflow on B335 
is to assume that the outflow has completely cleared out the dense material in a bipolar 
cone of constant opening angle. Figure 10 shows the radial behavior of the (H — K) colors 
from simulated observations we have made of an infalling model for four different values of 
the model outflow semi-opening angle a. (The appendix describes how the simulations are 
constructed.) Stars viewed through the model outflow cone are marked in red, and those 
not viewed through the outflow cone are marked in blue. The density distribution for 6 > a 
in each of the simulations is the infalling model that best fits the E(H — K) data beyond 
45° from the outflow axis (Fit B). As the model opening angle increases, the stars begin to 
separate into two distinct populations. Two effects are evident as the angle increases: (1) 
the population of detected stars with reduced colors becomes larger, and (2) the separation 
in the colors of the two populations grows larger. Manifest in the simulation with the largest 
opening angle (55°) is that the first effect is not only caused by the outflow cone occupying 
a larger fraction of the sky, but is also due to there being less column density for lines-of- 
sight passing through the cone, which allows a greater number of the background stars to 
be detected. Comparison of the radial behavior of the (H — K) colors in these simulations 
with the observations in Figure 7 shows that these simple models successfully reproduce the 
major features. The extent of the spread between the colors of the two populations of stars 
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in the NICMOS data seems to be reproduced best by the simulations with semi-opening 
angles of 35° and 45°. These angles also roughly reproduce the observed number of stars 
with significantly reduced colors. 

3.7. Including the Outflow (Fits D & E) 

The simple hollow cone model of the bipolar outflow allows the stars viewed closer to 
the outflow axis to be included in the x 2 analyses. The fact that the predicted (H — K) colors 
are a sensitive function of the model outflow opening angle means that we are able to fit for 
the opening angle of the flow. 

In both Fits D & E, we fit for the entire region within r < 100", with a model which has 
both the infall radius and outflow opening angle as free parameters, and with the color excess 
normalization for each model chosen to minimize the xl m the region beyond 45° from the 
outflow axis (as in Fit B). The only difference between Fit D and Fit E is that Fit E includes 
an extra 20% x E(H — K) dispersion term in &f lc to account for additional structure in the 
cloud, perhaps from the residual turbulent fluctuations. Table 1 lists the best-fit parameters 
for Fit D and Fit E, and they are nearly identical. In both Fits D & E, xl i s essentially an 
independent function of the infall radius and the outflow opening angle. Figure 11 shows 
a contour plot of the surface of xl f° r Fit D. The minimized value of xl is 4.02, with best 
fit values of the infall radius and outflow opening angle of R m f = 30 ± 3" (0.036 ± 0.004 pc) 
and a = 41 ± 4°, respectively. For Fit E, the contour plots and fitted parameters are 
essentially identical — infall radius i? in f = 26 ± 3" (0.031 ± 0.004 pc) and outflow opening 
angle a — 41 ± 2° — but the minimized xl is — 1 due to the extra dispersion. That the two 
fits agree so closely, particularly with regard to the infall radius, is interesting. The extra 
dispersion in Fit E reduces the weight of the most reddened (and most central) stars, which 
inevitably are those that constrain the fitting of the infall radius. 

The composite model with infall and bipolar outflow provides a consistent description 
of the density structure of B335. When the regions affected by the outflow are excluded 
from consideration, the inside-out collapse solution provided a good description of the data. 
Modeling the outflow as a hollowed-out bipolar cone allows inclusion of the stars from the 
outflow region in the x 2 analysis. Despite the obvious crudeness of this model for the outflow, 
the fits to the entire dataset (Fits D & E) constrain both the infall radius and outflow opening 
angle, and the fitted infall radii are entirely consistent with our earlier determination (Fit 
B) and with the values determined by molecular line studies. The fitted opening angle 
of the outflow is somewhat larger than indicated in studies of high velocity CO emission. 
Examination of the CO velocity structure (e.g. Hirano et al. 1988, 1992) suggests that the 
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flow is more collimated at higher velocities, and that the opening angles dervied from CO 
must be lower limits since the low velocity outflow becomes confused with the ambient globule 
emission (and is resolved out by interferometers). This effect may explain why the opening 
angle appears wider in extinction than in CO emission. Of course, the one parameter model 
fitted here simply may not provide an adequate description of the outflow geometry at this 
level of detail. 

Figure 12 presents a simulation of the models for Fit D and Fit E to afford another 
way to examine the success of these models, as well as the differences between them (the 
appendix describes the detailed recipe for constructing these simulations). The model for 
Fit D has an infall radius of 30" and an outflow opening angle of 41° (the normalization of 
the color excess is matched to the data by the x 2 analysis). This model reproduces the gross 
structure of the NICMOS data very successfully. However, the radial plot shows that the 
dispersion in the (H — K) colors at a given radius in the simulation is clearly significantly 
less than in the observations. The model for Fit E is nearly identical to Fit D except that 
the color excess has been multiplied by a Gaussian random variable with mean of unity and 
standard deviation 20%. The simulated data from the model for Fit E, with this additional 
dispersion to mimic the effects of residual turbulence, compares much more favorably with 
the observations in Figure 7. 



3.8. The Scaling Factor 

So far, we have been concerned only with the shape of the density profile and have 
ignored the free scaling factor applied to the color excess predicted by the various models 
during the fitting. For both the inside-out collapse models and the Bonnor-Ebert spheres, 
the normalization of the density distribution is entirely determined by the effective sound 
speed, measured to be 0.23 km s -1 (Zhou et al. 1990). This measurement, together with 
the distance to the cloud, gas-to-dust ratio, and near-infrared reddening law, make a unique 
prediction for the color excess normalization in each model. As listed in Table 1, a scaling 
factor T ~ 3-5 is required to match the models with the observed color excess. Figure 13 
shows a simulation of the spherically symmetric inside-out collapse model with the effective 
sound speed of 0.23 km s _1 , the standard distance of 250 pc (Tomita et al. 1979), and the 
infall radius of 25" derived by Choi et al. (1995). The simulation demonstrates the extent 
to which the model underestimates the observed reddening and extinction in B335. If this 
model were correct, then the NICMOS observations would have been sufficiently deep to see 
right through the center of the cloud, and the number of stars detected would have been 
far greater than observed. This therefore leads to the question: if the inside-out collapse 
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model or the Bonnor-Ebert sphere model is correct, then what explains the need for this 
large scaling factor between the predicted extinction and the measured values? 

The first possibility to consider is that the effective sound speed in B335 has been un- 
derestimated. To account for the T ~ 3-5 by a larger kinetic temperature in the initial 
state seems implausible — a temperature of T ~ 40-60 K would be needed. Multi-transition 
observations of NH3 show that the bulk of the dense gas in B335 has a kinetic temperature 
of 10-12 K (Menten et al. 1984), consistent with a balance between heating and cooling pro- 
cesses. To account for the factor with non-thermal motions also seems implausible, given the 
measured widths of molecular lines. A non-thermal contribution of Vturb — 0.35-0.45 km s" 1 
would be required, which is clearly at odds with a multitude of molecular line data. For 
example, Zhou et al. (1990) used observations of H 2 CO to determine a small turbulent con- 
tribution of fturb = 0.085 km s" 1 to the thermal sound speed. Similar results from a variety 
of different molecules are documented by Menten et al. (1984) and Frerking, Langer & Wil- 
son (1987). One might envision a form of non-thermal support that is not reflected in the 
observed molecular line-widths, perhaps because of a special viewing geometry. If magnetic 
support were primarily from a static field and not waves, and this field were oriented so that 
the motions were not observable along the line-of-sight, then it would be effectively unde- 
tectable. Although such a special situation seems contrived and unlikely, the orientation of 
B335 is special in that the outflow lies close to the plane of the sky, which may make such 
a configuration more plausible. 

Another possibility that affects the overall scaling lies in the outer boundary, i.e. where 
the cloud ends. Zhou et al. (1990) calculated the mass of B335 using an outer radius of 
0.15 pc based on their radio map and the 250 pc distance estimate. The molecular line 
profiles of dense gas tracers are not particularly sensitive to this assumption. We have 
adopted this value for the outer radius throughout the extinction modeling. Indeed, the 
column densities are also insensitive to this assumption since the density falls off steeply 
with radius. For a star at the infall radius of the best-fitting inside-out collapse model (Fit 
B), doubling the outer radius of the cloud increases the predicted extinction by only 7.5%, 
which is an unimportant amount with respect to the derived scaling factors. 

An overestimate of the distance to B335 would also result in predicted column densities 
that are too small. This distance uncertainty may be the dominant source of systematic 
error in the study. The distance of 250 pc was adopted from optical measurements of stellar 
reddening versus distance by Tomita et al. (1979), and this distance is uncertain at the 50% 
level. If B335 were closer, then the extinction would be larger at the same angular radius. To 
account for the full scaling factor would require that the distance be ~ 3-5 times closer, only 
about 50-80 pc, which seems unreasonable. Also, for a given level of observed extinction, 
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the inferred mass of the globule is proportional to the square of its distance (since the angle 
subtended by B335 is fixed). For the outer boundary assumption of Zhou et al. (1990), the 
mass inferred from our extinction data is 



The 14 M Q for a distance of 250 pc may be large for a small globule, but reducing the 
distance by the whole scaling factor lowers the mass to less than 1 M for the inside-out 
collapse model, and a little over 1 M for the Bonnor-Ebert sphere. This seems small given 
that the globule must be forming a star with a substantial fraction of a solar mass, and 
losing a lot of mass through outflow in the process. Furthermore, in the context of the 
inside-out collapse model, a reduction in the distance also results in a smaller mass for the 
protostellar core, as well as a shorter time since the onset of collapse (both are proportional 
to the distance for a given sound speed and infall radius). In short, a closer distance for 
B335 would ameliorate the discrepancy between the observed and predicted extinction, but 
it is not clear if this could account for the entire difference. 

A final uncertainty is the conversion from column density to color excess using the gas- 
to-dust ratio and reddening law. Significant variations in the reddening law have not been 
observed, and they are unlikely to explain much of the observed discrepancy in the color 
excess (although we briefly return to this point in the appendix). By contrast, the gas-to- 
dust ratio has been observed to vary in the Galaxy, with variations of ~ 2 not uncommon 
along different lines of sight (Bohlin et al. 1978). If the gas-to-dust ratio in B335 were 
lower than the standard value for the interstellar medium, or the reddening law steeper, 
then this would cause us to underestimate the extinction. Potentially, the full scaling factor 
could be accounted for by this uncertainty, or perhaps by some combination of a closer 
distance and a modified gas-to-dust ratio. For example, if the true distance to B335 were 
actually 100 pc, and the gas-to-dust ratio were about one half of the standard value (i.e. 
R — 1 x 10 21 cm -2 /mag), then these two factors would fully account for the discrepancy 
between observed extinction and the predictions of either model, and also give a reasonable 
mass for the globule of M ~ 2 M Q . If molecular hydrogen were to deplete onto dust grains, 
this could cause such a decrease in the gas-to-dust ratio at high density. However, H 2 is 
not thought to undergo significant freeze-out in this way (Sandford & Allamandola 1993), 
and a decrease in the gas-to-dust ratio is at odds with evidence that suggests an increase 
in high extinction regions (e.g. Bohlin et al. 1978, de Geus & Burton 1991, and Ciolek & 
Mouschovias 1996). 




(12) 
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3.9. Comparison with Dust Emission Studies 

The density distribution in protostellar cores may also be studied through dust emission. 
If the emission is observed at long enough wavelength that the Rayleigh- Jeans approximation 
applies, then, in the optically thin regime, the intensity at impact parameter b is proportional 
to the integral along the line-of-sight of the product of the density and the temperature. 
Assuming that the intensity, temperature and density obey power law distributions in radius, 
then the power law index of the density distribution p oc r _p is simply given by p = m + 1 — q, 
where m is the observed power law index of the intensity profile, and q the power law index of 
the temperature profile. The mass may also be derived from the radial intensity distribution, 
given an assumed mass opacity. Unfortunately, the mass opacity of dust at millimeter and 
submillimeter wavelengths is uncertain, by a factor of five or more (Ossenkopf & Henning 
1994). 

For B335, dust emission at submillimeter and millimeter wavelengths has been studied 
recently by Shirley et al. (2000) and Motte & Andre (2001). Shirley et al. (2000) constructed 
images at 850 /im and 450 /im at resolutions of 15" and 8", respectively, using SCUBA at the 
JCMT. At both wavelengths, these observations are sensitive to emission out to an angular 
radius of roughly 52". They adopt a temperature profile with power law index q = 0.4 
for a region where the envelope is optically thin to the bulk of the infrared radiation (e.g. 
Emerson 1988; Butner et al. 1990). They fit the circularly averaged maps at each wavelength, 
extending from inner cutoffs at angular radii of 24" for the 850 /im profile and 12" for the 
450 /im profile. They derive power law indices of p = 1.74 ±0.4 and p = 1.65 ±0.17 from the 
850 /im and 450 /im emission, respectively. They estimate a mass for B335 of 1.2 M using 
the total 850 /im flux in a 120" area and assuming a sphere of constant density. Motte and 
Andre (2001) constructed a 1.3 mm continuum image of B335 with 11" resolution using the 
IRAM 30m telescope, detecting emission out to an angular radius of 120". They adopt an 
isothermal temperature distribution and calculate a density power law index of p = 2.2 ±0.4. 
In this model, the total mass is 2.5 M . They note an asymmetry in the emission, with the 
emission being elongated perpendicular to the outflow axis. As in this extinction study, they 
perform fits to the off-outflow regions, and as a result increase the uncertainty in their fitted 
density index by 0.2 (included in the quoted uncertainty). 

Taking account the uncertainties, the density power law indices derived from the mil- 
limeter and submillimeter emission agree both with each other and with the results from 
dust extinction analysis (Fit A, p = 1.91 ±0.07). However, the extinction method provides a 
direct measurement of the distribution that suffers none of the added uncertainties associated 
with the temperature profile in the core. 
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4. Summary 

We present a near-infrared extinction study of the protostellar collapse candidate B335 
using deep integrations from HST/NICMOS supplemented by data from the W.M. Keck 
Observatory. In summary: 



1. We determine transformations between the NICMOS filters F160W and F222M and 
the standard H and K bands in the CIT photometric system for a range of (H — K) 
color from to 3. At the accuracy of our data, the derived transformations do not 
have significant slopes or color terms. 

2. The NICMOS mosaic image of the central 72" region of B335 shows a dramatic fall off 
in the number of stars detected towards the protostar location, where the extinction 
increases due to the central concentration of dense core material. The NICMOS mosaic 
image contains a total of 119 stars detected in both filters, the innermost of which is 
at a radius of 14". The photometry shows a steep gradient in the (H — K) colors 
toward the central protostar. In addition, there is a strong asymmetry in the colors 
superimposed on the overall gradient that coincides with the bipolar outflow. 

3. We compare a series of models of dense core stucture to the extinction data, includ- 
ing the previously proposed models of inside-out collapse derived from molecular line 
observations. The shape of the E(H — K) radial profile is well matched by the inside- 
out collapse model, in particular the density fall-off as r~ 2 for the envelope and a 
flattening within the purported infall zone. An unstable Bonnor-Ebert sphere with 
dimensionless outer radius £ max = 12.5 ± 2.6 provides an equally good description of 
the density profile, and is indistinguishable from the collapse model over the range 
in radius where there are stars to fit. The inside-out collapse model is not unique. 
The power law index for the density distribution derived from extinction is consis- 
tent with recent determinations from dust emission measurements at millimeter and 
submillimeter wavelengths. 

4. The bipolar outflow appears to be responsible for the dominant asymmetry in the 
extinction data. Stars whose line of sight fall within 45° of the outflow axis show lower 
extinction than those further away from the outflow. The effect of the outflow on the 
extinction data is reproduced very well by models that consider a hollowed out bipolar 
cone of constant opening angle. In the context of the inside-out collapse models, the 
observed asymmetry is too large to be explained by either rotation, or a weak magnetic 
field. 
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5. The entire extinction dataset can be well fitted by a model that includes infall, bipolar 
outflow, and an additional 20% x E(H — K) dispersion that mimics residual turbulent 
motions. The best fit outflow semi-opening angle is a = 41 ± 2° and the best fit 
infall radius is i?i n f = 26 ± 3". The fitted value of the infall radius is consistent with 
those derived from molecular line studies of B335, and supports the inside-out collapse 
interpretation of the density structure. The fitted opening angle for the bipolar outflow 
is somewhat larger than has been observed in high velocity CO emission, perhaps 
because the molecular line studies do not recover the full outflow width due to confusion 
with ambient emission at low velocities. 

6. The normalization of the observed color excess is a factor of T ~ 5 higher than the 
value predicted using the nominal 250 pc distance to B335, and the standard gas- 
to-dust ratio and near-infrared reddening law. The discepancy can conceivably be 
explained by a combination of a closer distance and a lower gas-to-dust ratio. 
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A. Simulated Extinction Observations 

This appendix contains the tools and general recipe for constructing the simulated 
observations of B335 presented in Sections 3.7 & 3.8. 

A.l. The Background Luminosity Function 

We construct a composite luminosity function for the background population by com- 
bining the data from fields away from the extinction of the B335 dense core. In deriving the 
luminosity function, the photometry for each field was cut at a very high completeness level, 
where the number of stars per unit magnitude range peaks. (This is a very conservative 
estimate of limiting magnitude.) To avoid any issues concerning how the stars should be 
binned, the cutoff was determined by finding where the gradient of the cumulative luminosity 



25 



function was steepest. For the dataset from the fifth NICMOS strip field, which includes 202 
detected stars at H and 85 detected stars at K, the cut-off magnitudes are 20.7 at H, and 
18.7 at K. For the combined dataset of the four Keck off-fields, which includes 242 detected 
stars at H and 246 detected stars at K, the cut-off magnitudes were 18.8 for both filters. 

The empirical H-band luminosity function was calculated in three parts out to the 
NICMOS high completeness limit (H = 20.7). For the bright end (H < 12.85), only stars from 
the NICMOS fifth strip field are included, as stars brighter than this limit were saturated in 
the Keck images. For the mid-range (12.85 < H < 18.8), stars from both datasets contribute 
to the function. For the faint end (18.8 < H < 20.7), above the Keck high completeness limit, 
only stars in the NICMOS fifth strip field were included. In each part, the number of stars 
was divided by the relevant observed area to build up the cumulative luminosity function. 
In the part where the NICMOS and Keck data overlap (12.85 < H < 18.8), the two datasets 
show good agreement. The K band luminosity function was constructed following the same 
procedure, and it is effectively the same as for the H band, offset by the mean color of 
H — K = 0.13 magnitudes. 



The central regions of B335 were observed to a deeper limiting magnitude than the 
strip or the 10' off fields, by about 4 magnitudes (H band) at the mosaic center, and about 
3 magnitudes at the edge (see Figure 2). Although the stars in these central regions of B335 
are highly extincted, the empirical background luminosity function must be extrapolated in 
order to simulate the deepest observations. Since the least extincted star in the NICMOS 
mosaic has H — K = 0.47 ± 0.08, corresponding to an extinction of = 1.3 ± 0.2, the 
luminosity function must be extrapolated by about 2 magnitudes at H in the simulations. 

The form of the H cumulative luminosity function suggests that a broken power law 
provides a good description. The function was fitted by minimizing the x 2 °f the fit to the 
bright-end of the luminosity function. The resulting function is given by: 



A. 1.1. Luminosity Function Extrapolation 



N(m < H) = 5.4 stars 




for H < 14.17 



for H > 14.17 



(Al) 
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A.2. A Recipe for Constucting Simulations 

Simulated observations of B335 were constructed by superimposing a model cloud on 
a field of background stars with the appropriate characteristics. The column density was 
calculated at the position of each background star and converted into two observable quan- 
tities: (1) an extinction at H band, and (2) a color excess E(H — K). These conversions are 
based on the gas-to-dust ratio, and reddening law. (A brief discussion on the uncertainties 
of these quantities is given in Section 2.1.) 

To summarize, the recipe for constructing simulated observations has the following steps: 



1. Generate a background star field by selecting randomly from the broken-power-law 
cumulative luminosity function in concert with the Gaussian color distribution of the 
background population. The broken-power-law luminosity function is extended to 
H = 25, just beyond the limiting magnitude at the center of the NICMOS mosaic (Note 
that model background stars intrinsically fainter than H = 23 are not recovered in the 
simulated observations of the fitted models from Table 1 — the function is extended to 
H = 25 simply to ensure completeness). Positions for each star are assigned randomly 
within a region 8' x 2', to encompass the area of the NICMOS observations. 

2. Select a model density distribution to describe B335. 

3. Calculate the column density for each star in the region. 

4. Convert the column density to an extinction at H band, An, and a color excess E(H— K) 
using an assumed gas-to-dust ratio, R, and a reddening law, An/Ay. The "standard" 
prescription is (Bohlin et al. 1978, Rieke & Lebofsky 1985): 

A 088 ( + H 2 ) \ ^xlCPcm^mag^ ( A H /A y \ 

Au = °- 88 UxltPcm-'A R ) \0A75- ) ' (A2) 

E(H-K) = A u / 2.78. (A3) 

(In the simulations, the reddening law was allowed to vary up to 15-25% steeper 
than the standard law in order to more exactly reproduce the observed starcounts. A 
fluctuation upward in the star counts at the one standard deviation level explains the 
inferred steeper reddening law.) 

5. Modify the magnitudes and colors of each background star appropriately. 



6. 



Simulate the NICMOS observations by applying the spatially dependent detection 
limits (Figure 2). 
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Table I. Summary of the x 2 Analyses Data 



Fit 


Model Density Distribution 


Fit Region 


Stars 


Fitted Model Parameter (s) 


xl 


A 


Scaled Power Law: p oc r~ p 


r < 100" ; 6 > 45° 


146 


p = 1.91 ±0.07 


3.37 


B 


Scaled Infall: p = •7 r pshu( r > Rini) 


r < 100" ; 9 > 45° 


146 


-R inf = 28 ± 3" ; T = 5.25 ± 0.13 


3.10 


C 


Scaled Bonnor-Ebcrt (R = 125") 


r < 100" ; 6 > 45° 


146 


£ max = 12.5 ± 2.6 ; T = 3.42 ± 0.10 


2.93 


D 


Scaled Infall, Hollow-Cone Outflow 


r < 100" 


209 


Rinf = 30 ± 3" ; a = 41 ± 4° ; T = 5.25 ± 0.13 


4.02 


E 


As Fit D with 20% Dispersion 


r < 100" 


209 


R inf = 26 ± 3" ; a = 41 ± 2° ; T = 4.77 ± 0.14 


1.04 



Note. — For each Fit, we consider only stars within 100" of the protostar location, to avoid noise-related issues near the edge of 
the cloud. Fits A-C consider only the regions that are well removed from the outflow, while Fits D & E are for the entire dataset 
within 100". Fit A is for a general power law density distribution. The remaining Fits use the Zhou et al. (1990) effective sound speed 
(a = 0.23 km s _1 ), with a free scaling applied to the predicted color excess, representing the many potential sources of systematic 
error in the conversion from color excess to density. Fit B uses the Shu inside-out collapse model, and Fit C uses a Bonnor-Ebert 
sphere. Fit D uses a combination of an infalling model with a model outflow that empties a cone of constant opening angle. Fit E 
uses the same model, but includes an extra dispersion in E(H — K) to account for extra small scale structure in the cloud. 
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Fig. 1. — Wide field optical view of the B335 region from the Digital Sky Survey. Superim- 
posed is a diagram of the observed fields from both NICMOS and Keck. 




Fig. 2. — Diagram showing the spatially dependent limiting magnitudes (in the CIT system) 
for the NICMOS observations. The quantity N indicates the number of exposures taken of 
that region. 
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Fig. 3.— Plot of NICMOS (F160W - F222M) color against Keck (H - K) color. Stars for 
which the Keck H or K measurement are of lower quality are marked as open circles. Also 
plotted is the derived transformation: (F160W - F222M) = (H - K) + 0.017. 
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Fig. 4. — Histogram of the (H — K) colors of the background population, with bins of width 
0.1 magnitudes. Also plotted is a Gaussian distribution with the same mean (H — K = 0.13) 
and standard deviation (cr(H — K) = 0.16), and normalized so that the area under the curves 
are equal. Stars observed at Keck with lower quality photometry as discussed in the text 
are not included in this histogram. 
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Fig. 5.— NICMOS F160W (left) and F222M (right) mosaic images of the central B335 
region. 
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Fig. 6.— NICMOS F160W (left) and F222M (right) images of the strip leading off B335. 
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Fig. 7. — Left: Pseudo-image of the NICMOS mosaic: the axes are spatial, magnitude 
determines a star's "size" , and (H — K) determines its color. The axes are in the frame of 
the NICMOS image. The mosaic contains 119 stars detected in both filters. The dashed line 
distinguishes the regions that are more than 45° from the outflow axis (the outflow axis has 
position angle roughly 100° in the pseudo-image). Right: (H — K) color against radius out 
to 180" from the center of B335. There is a steep gradient towards the protostar location. 
The symbol color indicates whether the star lies within 45° of the outflow axis. The outer 
edge of B335 occurs at a radius of approximately 125". 




Fig. 8. — Schematic diagram showing the scale and coverage of the NICMOS observations 
relative to important scales of B335. The region beyond 45° from the outflow axis is marked. 
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Fig. 9. — (H — K) vs. radius for the NICMOS stars beyond 45° from the outflow axis. 
Overplotted are predicted curves for the best fitting infall model (Fit B) and the unstable 
Bonnor-Ebert sphere (Fit C), assuming a constant background color of 0.13. The two curves 
are indistinguishable over the range in radius where there are stars to fit. 
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Fig. 10. — (H — K) vs. radius for simulations of an infalling model with outflow cones of 
various outflow semi-opening angles, a. 
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Fig. 11. — Contours of the surface of xt f° r models with both infall and outflow, where the 
outflow has been assumed to clear out a cone of opening angle (a). The minimum in the 
reduced \ 2 is 4.02 and corresponds to fitted parameters of R[ n f = 30 ± 3" (0.034 ± 0.04 pc at 
a distance of 250 pc), and a = 41 ± 4°. The position of the fitted paramters is marked with 
a diamond. The best fitting model is labelled Fit D in Table 1. 



-42- 





Fig. 12. — A simulation of the models of Fit D (Top) and Fit E (Bottom). For each model: 
(Left) Pseudo-image of the simulated NICMOS mosaic image; (Right) (H — K) color vs. 
radius out to 180" from the center. 
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Fig. 13. — A simulation of the spherically symmetric inside-out collapse model with infall 
radius 25" (Choi et al. 1995), and normalization determined by the sound speed of a = 
0.23 km s _1 (Zhou et al. 1990) and standard distance of 250 pc (Tomita et al. 1979). Left: 
Pseudo-image of the simulated NICMOS mosaic image. Right: (H — K) color vs. radius out 
to 180" from the center. 



